
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The impact of divorce laws on the equilibrium in the marriage market.
% Ana Reynoso
% April 2024
%
% This file inputs the cleaned PSID data and creates moments used in estimation 
% and their confidence intervals.
%
% Data: PSID 1968-1992
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

*%------------------------ Replication Path -------------------------------%

clear all
clear
set mem 1000M
set maxvar 15000
set more off, permanently

*%--- Indicate location of Replication folder:
*global replication_location "C:\update_with_your_path"

*%------------------------ Preliminaries ----------------------------------%

cd "$replication_location\Data\Inputs"

*%------------------------ Data moments ----------------------------------%
qui{
*** Market 1: Northeast
clear all
set more off, permanently

use data_for_moments, clear
keep if USCBregion==1

gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

*- Matching patterns

*-- Number of couples by couple type
gen nbr_couple_type=.
forvalues g=1(1)9{
qui sum couple_type if couple_type==`g'
replace nbr_couple_type=r(N) if couple_type==`g'
}

*-- Number of singles by educ
forvalues s=10(1)15{
qui sum couple_type if couple_type==`s'
replace nbr_couple_type=r(N) if couple_type==`s'
}

*-- Number of females
gen nbr_fem=.
qui sum theres_female if theres_female==1
replace nbr_fem=r(N)

*-- Number of females by educ type
bysort female_type: egen nbr_females_bytype=sum(theres_female)

*-- Number of males by educ type
bysort male_type: egen nbr_males_bytype=sum(theres_male)

*-- Relative measures of males and females
gen hs_fem=nbr_females_bytype/nbr_fem if female_type==1
gen sc_fem=nbr_females_bytype/nbr_fem if female_type==2
gen cp_fem=nbr_females_bytype/nbr_fem if female_type==3

gen hs_male=nbr_males_bytype/nbr_fem if male_type==1
gen sc_male=nbr_males_bytype/nbr_fem if male_type==2
gen cp_male=nbr_males_bytype/nbr_fem if male_type==3

*-- Frequency of couples relative to number of households with females
gen mu_f=.
forvalues g=1(1)9{
replace mu_f= nbr_couple_type/nbr_fem if couple_type==`g'
}
*-- Frequency of singles by education relative to number of people of that education
gen prob_s=.
forvalues s=10(1)12{
replace prob_s= nbr_couple_type/nbr_females_bytype if couple_type==`s'
}
forvalues s=13(1)15{
replace prob_s= nbr_couple_type/nbr_males_bytype if couple_type==`s'
}

*- stay-at-home wife pooled over the life cycle
gen hw_lc=.
forvalues g=1(1)9{
qui sum sahw if couple_type==`g'
replace hw_lc=r(mean) if couple_type==`g'
}

*- divorce probability
gen div_hazard=.
forvalues g=1(1)9{
qui sum div_hazard_dummy if couple_type==`g'
replace div_hazard=r(mean) if couple_type==`g'
}

*-save
keep couple_type div_hazard hw_lc mu_f prob_s hs_fem- cp_male 
duplicates drop
sort couple_type
*-- Measures of number of people per type relative to number of females
gen frac_fem_type=hs_fem
replace frac_fem_type=sc_fem if couple_type>=4
replace frac_fem_type=cp_fem if couple_type>=7
gen frac_male_type=hs_male
replace frac_male_type=sc_male if sc_male!=.
replace frac_male_type=cp_male if cp_male!=.
rename frac_fem_type frac_fem_type_mkt1
rename frac_male_type frac_male_type_mkt1

save data_frac_mkt1.dta, replace

keep couple_type mu_f prob_s hw_lc div_h
gen mu_f2=mu_f
stack mu_f mu_f2 prob_s hw_lc div_h, into(moments)
drop if moments==.
drop _stack
gen mom_nbr=_n
rename moments moments_mkt1
sort mom_nbr

save data_moments_mkt1.dta, replace

*** Market 2: West and Midwest

clear all
set more off, permanently

use data_for_moments, clear
keep if USCBregion==2

gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

*- Matching patterns
	
*-- Number of couples by couple type
gen nbr_couple_type=.
forvalues g=1(1)9{
qui sum couple_type if couple_type==`g'
replace nbr_couple_type=r(N) if couple_type==`g'
}

*-- Number of singles by educ
forvalues s=10(1)15{
qui sum couple_type if couple_type==`s'
replace nbr_couple_type=r(N) if couple_type==`s'
}

*-- Number of females
gen nbr_fem=.
qui sum theres_female if theres_female==1
replace nbr_fem=r(N)

*-- Number of females by educ type
bysort female_type: egen nbr_females_bytype=sum(theres_female)

*-- Number of males by educ type
bysort male_type: egen nbr_males_bytype=sum(theres_male)

*-- Relative measures of males and females (with respect to measure of females)
gen hs_fem=nbr_females_bytype/nbr_fem if female_type==1
gen sc_fem=nbr_females_bytype/nbr_fem if female_type==2
gen cp_fem=nbr_females_bytype/nbr_fem if female_type==3
gen hs_male=nbr_males_bytype/nbr_fem if male_type==1
gen sc_male=nbr_males_bytype/nbr_fem if male_type==2
gen cp_male=nbr_males_bytype/nbr_fem if male_type==3

*-- Frequency of couples relative to number of households with females
gen mu_f=.
forvalues g=1(1)9{
replace mu_f= nbr_couple_type/nbr_fem if couple_type==`g'
}

*-- Frequency of singles by education relative to number of people of that education
gen prob_s=.
forvalues s=10(1)12{
replace prob_s= nbr_couple_type/nbr_females_bytype if couple_type==`s'
}
forvalues s=13(1)15{
replace prob_s= nbr_couple_type/nbr_males_bytype if couple_type==`s'
}

*- stay-at-home wife pooled over the life cycle
gen hw_lc=.
forvalues g=1(1)9{
qui sum sahw if couple_type==`g'
replace hw_lc=r(mean) if couple_type==`g'
}

*- divorce probability
gen div_hazard=.
forvalues g=1(1)9{
qui sum div_hazard_dummy if couple_type==`g'
replace div_hazard=r(mean) if couple_type==`g'
}

*- save
keep couple_type div_hazard hw_lc mu_f prob_s hs_fem- cp_male 
duplicates drop
sort couple_type

gen frac_fem_type=hs_fem
replace frac_fem_type=sc_fem if couple_type>=4
replace frac_fem_type=cp_fem if couple_type>=7
gen frac_male_type=hs_male
replace frac_male_type=sc_male if sc_male!=.
replace frac_male_type=cp_male if cp_male!=.
rename frac_fem_type frac_fem_type_mkt2
rename frac_male_type frac_male_type_mkt2

save data_frac_mkt2.dta, replace

keep couple_type mu_f prob_s hw_lc div_h
gen mu_f2=mu_f
stack mu_f mu_f2 prob_s hw_lc div_h, into(moments)
drop if moments==.
drop _stack
gen mom_nbr=_n
rename moments moments_mkt2
sort mom_nbr

save data_moments_mkt2.dta, replace

*** Market 3: South Atlantic 
clear all
set more off, permanently

use data_for_moments, clear
keep if USCBregion==3
gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

*- Matching patterns
	
*-- Number of couples by couple type
gen nbr_couple_type=.
forvalues g=1(1)9{
qui sum couple_type if couple_type==`g'
replace nbr_couple_type=r(N) if couple_type==`g'
}

*-- Number of singles by educ
forvalues s=10(1)15{
qui sum couple_type if couple_type==`s'
replace nbr_couple_type=r(N) if couple_type==`s'
}

*-- Number of females
gen nbr_fem=.
qui sum theres_female if theres_female==1
replace nbr_fem=r(N)

*-- Number of females by educ type
bysort female_type: egen nbr_females_bytype=sum(theres_female)

*-- Number of males by educ type
bysort male_type: egen nbr_males_bytype=sum(theres_male)

*-- Relative measures of males and females (with respect to measure of females)
gen hs_fem=nbr_females_bytype/nbr_fem if female_type==1
gen sc_fem=nbr_females_bytype/nbr_fem if female_type==2
gen cp_fem=nbr_females_bytype/nbr_fem if female_type==3
gen hs_male=nbr_males_bytype/nbr_fem if male_type==1
gen sc_male=nbr_males_bytype/nbr_fem if male_type==2
gen cp_male=nbr_males_bytype/nbr_fem if male_type==3

*-- Frequency of couples relative to number of households with females
gen mu_f=.
forvalues g=1(1)9{
replace mu_f= nbr_couple_type/nbr_fem if couple_type==`g'
}

*-- Frequency of singles by education relative to number of people of that education
gen prob_s=.
forvalues s=10(1)12{
replace prob_s= nbr_couple_type/nbr_females_bytype if couple_type==`s'
}
forvalues s=13(1)15{
replace prob_s= nbr_couple_type/nbr_males_bytype if couple_type==`s'
}

*- stay-at-home wife pooled over the life cycle
gen hw_lc=.
forvalues g=1(1)9{
qui sum sahw if couple_type==`g'
replace hw_lc=r(mean) if couple_type==`g'
}

*- divorce probability
gen div_hazard=.
forvalues g=1(1)9{
qui sum div_hazard_dummy if couple_type==`g'
replace div_hazard=r(mean) if couple_type==`g'
}

*-save

keep couple_type div_hazard hw_lc mu_f prob_s hs_fem- cp_male 
duplicates drop
sort couple_type

gen frac_fem_type=hs_fem
replace frac_fem_type=sc_fem if couple_type>=4
replace frac_fem_type=cp_fem if couple_type>=7
gen frac_male_type=hs_male
replace frac_male_type=sc_male if sc_male!=.
replace frac_male_type=cp_male if cp_male!=.
rename frac_fem_type frac_fem_type_mkt3
rename frac_male_type frac_male_type_mkt3

save data_frac_mkt3.dta, replace

keep couple_type mu_f prob_s hw_lc div_h
gen mu_f2=mu_f
stack mu_f mu_f2 prob_s hw_lc div_h, into(moments)
drop if moments==.
drop _stack
gen mom_nbr=_n
rename moments moments_mkt3
sort mom_nbr

save data_moments_mkt3.dta, replace

*** Market 4: South Central 
clear all
set more off, permanently

use data_for_moments, clear
keep if USCBregion==4
gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

*- Matching patterns:	
*-- Number of couples by couple type

gen nbr_couple_type=.
forvalues g=1(1)9{
qui sum couple_type if couple_type==`g'
replace nbr_couple_type=r(N) if couple_type==`g'
}

*-- Number of singles by educ
forvalues s=10(1)15{
qui sum couple_type if couple_type==`s'
replace nbr_couple_type=r(N) if couple_type==`s'
}

*-- Number of females
gen nbr_fem=.
qui sum theres_female if theres_female==1
replace nbr_fem=r(N)

*-- Number of females by educ type
bysort female_type: egen nbr_females_bytype=sum(theres_female)

*-- Number of males by educ type
bysort male_type: egen nbr_males_bytype=sum(theres_male)

*-- Relative measures of males and females (with respect to measure of females)
gen hs_fem=nbr_females_bytype/nbr_fem if female_type==1
gen sc_fem=nbr_females_bytype/nbr_fem if female_type==2
gen cp_fem=nbr_females_bytype/nbr_fem if female_type==3
gen hs_male=nbr_males_bytype/nbr_fem if male_type==1
gen sc_male=nbr_males_bytype/nbr_fem if male_type==2
gen cp_male=nbr_males_bytype/nbr_fem if male_type==3

*-- Frequency of couples relative to number of households with females
gen mu_f=.
forvalues g=1(1)9{
replace mu_f= nbr_couple_type/nbr_fem if couple_type==`g'
}

*-- Frequency of singles by education relative to number of people of that education
gen prob_s=.
forvalues s=10(1)12{
replace prob_s= nbr_couple_type/nbr_females_bytype if couple_type==`s'
}
forvalues s=13(1)15{
replace prob_s= nbr_couple_type/nbr_males_bytype if couple_type==`s'
}

*- stay-at-home wife pooled over the life cycle
gen hw_lc=.
forvalues g=1(1)9{
qui sum sahw if couple_type==`g'
replace hw_lc=r(mean) if couple_type==`g'
}

*- divorce probability
gen div_hazard=.
forvalues g=1(1)9{
qui sum div_hazard_dummy if couple_type==`g'
replace div_hazard=r(mean) if couple_type==`g'
}

*- save
keep couple_type div_hazard hw_lc mu_f prob_s hs_fem- cp_male 
duplicates drop
sort couple_type
gen frac_fem_type=hs_fem
replace frac_fem_type=sc_fem if couple_type>=4
replace frac_fem_type=cp_fem if couple_type>=7
gen frac_male_type=hs_male
replace frac_male_type=sc_male if sc_male!=.
replace frac_male_type=cp_male if cp_male!=.
rename frac_fem_type frac_fem_type_mkt4
rename frac_male_type frac_male_type_mkt4

save data_frac_mkt4.dta, replace

keep couple_type mu_f prob_s hw_lc div_h
gen mu_f2=mu_f
stack mu_f mu_f2 prob_s hw_lc div_h, into(moments)
drop if moments==.
drop _stack
gen mom_nbr=_n
rename moments moments_mkt4
sort mom_nbr

save data_moments_mkt4.dta, replace

*** All moments

use data_frac_mkt1, clear
keep couple_type frac*
merge couple_type using data_frac_mkt2
keep couple_type frac*
sort couple_type
merge couple_type using data_frac_mkt3
keep couple_type frac*
sort couple_type
merge couple_type using data_frac_mkt4
keep couple_type frac*

outsheet frac_fem_type* if couple_type<=9 using "data_f.csv" , comma replace 
outsheet frac_male_type* if couple_type<=9 using "data_m.csv" , comma replace

use data_moments_mkt1, clear
merge mom_nbr using data_moments_mkt2, sort
drop _merge
sort mom_nbr
merge mom_nbr using data_moments_mkt3, sort
drop _merge
sort mom_nbr
merge mom_nbr using data_moments_mkt4, sort
drop mom_nbr _merge

outsheet using "data_moments.csv", comma replace

}

*%------------------------ Bootstrapped data moments and CIs ------------------%

qui{
*** Market 1: Northeast

clear all
set more off, permanently

use data_for_moments, clear
keep if USCBregion == 1

gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

cap program drop prog_momentboot
program prog_momentboot, eclass
tempname ests
	mat `ests'=J(1,15+1+3+3+9+9,.)
*- Number of couples by couple type
forvalues g=1(1)15{
qui sum couple_type if couple_type==`g'
mat `ests'[1,`g']=r(N)
}
*- Number of females
qui sum theres_female if theres_female==1
mat `ests'[1,1+15]=r(N)
*- Number of females by type
forvalues i=1(1)3{
qui sum theres_female if theres_female==1 & female_type==`i'
mat `ests'[1,`i'+15+1]=r(N)
}
*- Number of males by type
forvalues i=1(1)3{
qui sum theres_male if theres_male==1 & male_type==`i'
mat `ests'[1,`i'+15+1+3]=r(N)
}
* Stay-at-home wife 
forvalues g=1/9{
qui sum sahw if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3]=r(mean)
}
* Divorce 
forvalues g=1/9{
qui sum div_hazard_dummy if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3+9]=r(mean)
}
	ereturn post `ests'
end

prog_momentboot
bootstrap _b, reps(1000) seed(1218) saving(momentboot_mat_1,replace): prog_momentboot

clear
use momentboot_mat_1
forvalues i=1/9{
gen mu_f_`i'=_b_c`i'/_b_c16
}
forvalues i=1/9{
gen mu_f2_`i'=mu_f_`i'
}
gen single_f_1=_b_c10/_b_c17
gen single_f_2=_b_c11/_b_c18
gen single_f_3=_b_c12/_b_c19
gen single_m_1=_b_c13/_b_c20
gen single_m_2=_b_c14/_b_c21
gen single_m_3=_b_c15/_b_c22
drop  _b_c1- _b_c22
order  mu_f_1 - single_m_3 
forvalues j=1/9{
rename mu_f_`j' mu_f_`j'_1
rename mu_f2_`j' mu_f2_`j'_1
}
forvalues j=1/3{
rename single_f_`j' single_f_`j'_1
rename single_m_`j' single_m_`j'_1
}
forvalues j=23/31{
rename _b_c`j' hw_`j'_1
}
forvalues j=32/40{
rename _b_c`j' dh_`j'_1
}
gen rep=_n
sort rep

save moments_boot_1.dta, replace

*** Market 2: West and Midwest

clear all
set more off, permanently

use data_for_moments, clear
keep if USCBregion == 2

gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

cap program drop prog_momentboot
program prog_momentboot, eclass

tempname ests
	mat `ests'=J(1,15+1+3+3+9+9,.)
*- Number of couples by couple type
forvalues g=1(1)15{
qui sum couple_type if couple_type==`g'
mat `ests'[1,`g']=r(N)
}
*- Number of females
qui sum theres_female if theres_female==1
mat `ests'[1,1+15]=r(N)
*- Number of females by type
forvalues i=1(1)3{
qui sum theres_female if theres_female==1 & female_type==`i'
mat `ests'[1,`i'+15+1]=r(N)
}
*- Number of males by type
forvalues i=1(1)3{
qui sum theres_male if theres_male==1 & male_type==`i'
mat `ests'[1,`i'+15+1+3]=r(N)
}
* Stay-at-home wife
forvalues g=1/9{
qui sum sahw if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3]=r(mean)
}
* Divorce 
forvalues g=1/9{
qui sum div_hazard_dummy if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3+9]=r(mean)
}
	ereturn post `ests'
end
prog_momentboot
bootstrap _b, reps(1000) seed(1218) saving(momentboot_mat_2,replace): prog_momentboot
clear
use momentboot_mat_2
forvalues i=1/9{
gen mu_f_`i'=_b_c`i'/_b_c16
}
forvalues i=1/9{
gen mu_f2_`i'=mu_f_`i'
}
gen single_f_1=_b_c10/_b_c17
gen single_f_2=_b_c11/_b_c18
gen single_f_3=_b_c12/_b_c19
gen single_m_1=_b_c13/_b_c20
gen single_m_2=_b_c14/_b_c21
gen single_m_3=_b_c15/_b_c22
drop  _b_c1- _b_c22
order  mu_f_1 - single_m_3 

forvalues j=1/9{
rename mu_f_`j' mu_f_`j'_2
rename mu_f2_`j' mu_f2_`j'_2
}
forvalues j=1/3{
rename single_f_`j' single_f_`j'_2
rename single_m_`j' single_m_`j'_2
}
forvalues j=23/31{
rename _b_c`j' hw_`j'_2
}
forvalues j=32/40{
rename _b_c`j' dh_`j'_2
}
gen rep=_n
sort rep

save moments_boot_2.dta, replace

*** Market 3: South Atlantic 
clear all
set more off, permanently
use data_for_moments, clear
keep if USCBregion == 3
gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

cap program drop prog_momentboot
program prog_momentboot, eclass

tempname ests
	mat `ests'=J(1,15+1+3+3+9+9,.)  
							  
*- Number of couples by couple type
forvalues g=1(1)15{
qui sum couple_type if couple_type==`g'
mat `ests'[1,`g']=r(N)
}
*- Number of females
qui sum theres_female if theres_female==1
mat `ests'[1,1+15]=r(N)
*- Number of females by type
forvalues i=1(1)3{
qui sum theres_female if theres_female==1 & female_type==`i'
mat `ests'[1,`i'+15+1]=r(N)
}
*- Number of males by type

forvalues i=1(1)3{
qui sum theres_male if theres_male==1 & male_type==`i'
mat `ests'[1,`i'+15+1+3]=r(N)
}
* Stay-at-home wife 
forvalues g=1/9{
qui sum sahw if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3]=r(mean)
}
* Divorce 
forvalues g=1/9{
qui sum div_hazard_dummy if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3+9]=r(mean)
}
	ereturn post `ests'
end

prog_momentboot
bootstrap _b, reps(1000) seed(1218) saving(momentboot_mat_3,replace): prog_momentboot
clear
use momentboot_mat_3

forvalues i=1/9{
gen mu_f_`i'=_b_c`i'/_b_c16
}
forvalues i=1/9{
gen mu_f2_`i'=mu_f_`i'
}
gen single_f_1=_b_c10/_b_c17
gen single_f_2=_b_c11/_b_c18
gen single_f_3=_b_c12/_b_c19
gen single_m_1=_b_c13/_b_c20
gen single_m_2=_b_c14/_b_c21
gen single_m_3=_b_c15/_b_c22
drop  _b_c1- _b_c22
order  mu_f_1 - single_m_3 
forvalues j=1/9{
rename mu_f_`j' mu_f_`j'_3
rename mu_f2_`j' mu_f2_`j'_3
}
forvalues j=1/3{
rename single_f_`j' single_f_`j'_3
rename single_m_`j' single_m_`j'_3
}
forvalues j=23/31{
rename _b_c`j' hw_`j'_3
}
forvalues j=32/40{
rename _b_c`j' dh_`j'_3
}
gen rep=_n
sort rep

save moments_boot_3.dta, replace

*** Market 4: South Central 
clear all
set more off, permanently
use data_for_moments, clear
keep if USCBregion == 4

gen div_hazard_dummy=ever_divorced if couple_marries==1
gen theres_male=1 if male_type!=.
gen theres_female=1 if female_type!=.

cap program drop prog_momentboot
program prog_momentboot, eclass

tempname ests
	mat `ests'=J(1,15+1+3+3+9+9,.) 
*- Number of couples by couple type
forvalues g=1(1)15{
qui sum couple_type if couple_type==`g'
mat `ests'[1,`g']=r(N)
}
*- Number of females
qui sum theres_female if theres_female==1
mat `ests'[1,1+15]=r(N)
*- Number of females by type
forvalues i=1(1)3{
qui sum theres_female if theres_female==1 & female_type==`i'
mat `ests'[1,`i'+15+1]=r(N)
}
*- Number of males by type
forvalues i=1(1)3{
qui sum theres_male if theres_male==1 & male_type==`i'
mat `ests'[1,`i'+15+1+3]=r(N)
}
* Stay-at-home wife 
forvalues g=1/9{
qui sum sahw if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3]=r(mean)
}
* Divorce 
forvalues g=1/9{
qui sum div_hazard_dummy if couple_type==`g'
mat `ests'[1,`g'+15+1+3+3+9]=r(mean)
}
	ereturn post `ests'
end

prog_momentboot
bootstrap _b, reps(1000) seed(1218) saving(momentboot_mat_4,replace): prog_momentboot 
clear
use momentboot_mat_4

forvalues i=1/9{
gen mu_f_`i'=_b_c`i'/_b_c16
}
forvalues i=1/9{
gen mu_f2_`i'=mu_f_`i'
}
gen single_f_1=_b_c10/_b_c17
gen single_f_2=_b_c11/_b_c18
gen single_f_3=_b_c12/_b_c19
gen single_m_1=_b_c13/_b_c20
gen single_m_2=_b_c14/_b_c21
gen single_m_3=_b_c15/_b_c22
drop  _b_c1- _b_c22
order  mu_f_1 - single_m_3 
forvalues j=1/9{
rename mu_f_`j' mu_f_`j'_4
rename mu_f2_`j' mu_f2_`j'_4
}
forvalues j=1/3{
rename single_f_`j' single_f_`j'_4
rename single_m_`j' single_m_`j'_4
}
forvalues j=23/31{
rename _b_c`j' hw_`j'_4
}
forvalues j=32/40{
rename _b_c`j' dh_`j'_4
}
gen rep=_n
sort rep

save moments_boot_4.dta, replace

*** Append all markets together

use moments_boot_1.dta, clear
merge rep using moments_boot_2
drop _merge
merge rep using moments_boot_3, sort
drop _merge
merge rep using moments_boot_4, sort
drop rep
drop _merge
outsheet using "moments_boot.csv", comma replace

*** Confidence intervals of the data 

use moments_boot_1.dta, clear
forvalues j=1/4{     
use moments_boot_`j'.dta, clear		
forvalues i=1/9{
qui sum mu_f_`i'_`j',d
gen MUF_l_`i'=r(p5)
replace MUF_l_`i'=round(MUF_l_`i',.0001)
gen MUF_u_`i'=r(p95)
replace MUF_u_`i'=round(MUF_u_`i',.0001)
}
forvalues i=1/3{
qui sum single_f_`i'_`j',d
gen Sf_l_`i'=r(p5)
replace Sf_l_`i'=round(Sf_l_`i',.0001)
gen Sf_u_`i'=r(p95)
replace Sf_u_`i'=round(Sf_u_`i',.0001)
}
forvalues i=1/3{
qui sum single_m_`i'_`j',d
gen Sm_l_`i'=r(p5)
replace Sm_l_`i'=round(Sm_l_`i',.0001)
gen Sm_u_`i'=r(p95)
replace Sm_u_`i'=round(Sm_u_`i',.0001)
}
forvalues i=23/31{
qui sum hw_`i'_`j',d
gen hw_l_`i'=r(p5)
replace hw_l_`i'=round(hw_l_`i',.0001)
gen hw_u_`i'=r(p95)
replace hw_u_`i'=round(hw_u_`i',.0001)
}
forvalues i=32/40{
qui sum dh_`i'_`j',d
gen d_l_`i'=r(p5)
replace d_l_`i'=round(d_l_`i',.0001)
gen d_u_`i'=r(p95)
replace d_u_`i'=round(d_u_`i',.0001)
}
drop mu_f_1_`j'-dh_40_`j'
keep if rep==1
save moments_ci_`j'.dta, replace
}     
use moments_ci_1.dta, clear
append using moments_ci_2
append using moments_ci_3
append using moments_ci_4
drop rep

outsheet using "moments_boot_ci.csv", comma replace

}

*%---------- Earnings regressions bootstrapped data moments and CIs------------%

qui{
	
clear all
set more off, permanently
u data_for_earnings-moments.dta
sort couple hh_id sample_person year

*- Convert income into real values (from Voena (2015) code)
foreach var in ly wly {
	gen real_`var'=.
	replace real_`var'=`var'*3.739003 if year==1968
	replace real_`var'=`var'*3.571429 if year==1969
	replace real_`var'=`var'*3.364116 if year==1970
	replace real_`var'=`var'*3.195488 if year==1971
	replace real_`var'=`var'*3.09466 if year==1972
	replace real_`var'=`var'*2.985949 if year==1973
	replace real_`var'=`var'*2.724359 if year==1974
	replace real_`var'=`var'*2.437859 if year==1975
	replace real_`var'=`var'*2.284946 if year==1976
	replace real_`var'=`var'*2.172061 if year==1977
	replace real_`var'=`var'*2.033493 if year==1978
	replace real_`var'=`var'*1.861314 if year==1979
	replace real_`var'=`var'*1.634615 if year==1980
	replace real_`var'=`var'*1.462156 if year==1981
	replace real_`var'=`var'*1.350636 if year==1982
	replace real_`var'=`var'*1.302349 if year==1983
	replace real_`var'=`var'*1.248776 if year==1984
	replace real_`var'=`var'*1.206244 if year==1985
	replace real_`var'=`var'* 1.2560146 if year==1986
	replace real_`var'=`var'* 1.2544524 if year==1987
	replace real_`var'=`var'*1.099138 if year==1988
	replace real_`var'=`var'*1.05198 if year==1989
	replace real_`var'=`var'*1 if year==1990
	replace real_`var'=`var'*0.9465479 if year==1991
	replace real_`var'=`var'*0.9219089 if year==1992
	replace real_`var'=`var'*0.8928571 if year==1993
	drop `var'
	rename real_`var' `var'
	}

*- Create variables
gen hrw=ly/hours 
gen lhrw=log(hrw)
gen l_ly=log(ly)
gen fem_ly=.
replace fem_ly=ly if gender==2
replace fem_ly=wly if s_gender==2
gen lfem_ly=log(fem_ly)

gen hw2=.
replace hw2=0 if lfem_ly!=.
replace hw2=1 if lfem_ly==. &(gender==2 | s_gender==2)
replace hw2=1 if s_hrs_wk_annual<520 & s_gender==2
replace hw2=1 if hrs_wk_annual<520 & gender==2

gen female_works=.
replace female_works=1-hw2 if hw2!=.
bysort woman (year) : gen X= sum(female_works)-female_works 
gen sqr_X= X*X

gen hwft=.
replace hwft=0 if lfem_ly!=.
replace hwft=1 if lfem_ly==. &(gender==2 | s_gender==2)
bysort woman (year) : gen K= sum(hwft)-hwft 
replace K=0 if couple==0 & man!=.

*- Save 
gen sqrage = maleage^2
gen sqrfage = femage^2
save data_for_earnings-moments-allmkts.dta, replace

*- Earnings moments by market

*---Mkt1
u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==1
forvalues i=1(1)3{
reg l_ly K maleage sqrage if man!=.  & maleduc_cat==`i' & couple==1, robust
gen bm0_`i'=_b[_cons]
gen bm1_`i'=_b[maleage]
gen bm2_`i'=_b[sqrage]
gen bm3_`i'=_b[K]
}
forvalues i=1(1)3{
reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==`i' & couple==1, robust
gen bf0_`i'=_b[_cons]
gen bf1_`i'=_b[X]
gen bf2_`i'=_b[sqr_X]
}

keep bm* bf*
duplicates drop
stack bm0_1-bf2_3, into(moments_mkt1)
drop _stack
gen mom_nbr=_n
sort mom_nbr

save reg_moments_mkt1.dta, replace

*---Mkt2
u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==2
forvalues i=1(1)3{
reg l_ly K maleage sqrage if man!=.  & maleduc_cat==`i' & couple==1, robust
gen bm0_`i'=_b[_cons]
gen bm1_`i'=_b[maleage]
gen bm2_`i'=_b[sqrage]
gen bm3_`i'=_b[K]
}
forvalues i=1(1)3{
reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==`i' & couple==1, robust
gen bf0_`i'=_b[_cons]
gen bf1_`i'=_b[X]
gen bf2_`i'=_b[sqr_X]
}

keep bm* bf*
duplicates drop
stack bm0_1-bf2_3, into(moments_mkt2)
drop _stack
gen mom_nbr=_n
sort mom_nbr

save reg_moments_mkt2.dta, replace

*---Mkt3
u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==3
forvalues i=1(1)3{
reg l_ly K maleage sqrage if man!=.  & maleduc_cat==`i' & couple==1, robust
gen bm0_`i'=_b[_cons]
gen bm1_`i'=_b[maleage]
gen bm2_`i'=_b[sqrage]
gen bm3_`i'=_b[K]
}
forvalues i=1(1)3{
reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==`i' & couple==1, robust
gen bf0_`i'=_b[_cons]
gen bf1_`i'=_b[X]
gen bf2_`i'=_b[sqr_X]
}

keep bm* bf*
duplicates drop
stack bm0_1-bf2_3, into(moments_mkt3)
drop _stack
gen mom_nbr=_n
sort mom_nbr

save reg_moments_mkt3.dta, replace

*---Mkt4
u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==4
forvalues i=1(1)3{
reg l_ly K maleage sqrage if man!=.  & maleduc_cat==`i' & couple==1, robust
gen bm0_`i'=_b[_cons]
gen bm1_`i'=_b[maleage]
gen bm2_`i'=_b[sqrage]
gen bm3_`i'=_b[K]
}
forvalues i=1(1)3{
reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==`i' & couple==1, robust
gen bf0_`i'=_b[_cons]
gen bf1_`i'=_b[X]
gen bf2_`i'=_b[sqr_X]
}

keep bm* bf*
duplicates drop
stack bm0_1-bf2_3, into(moments_mkt4)
drop _stack
gen mom_nbr=_n
sort mom_nbr

save reg_moments_mkt4.dta, replace

*---Append all markets together

use reg_moments_mkt1, clear
merge mom_nbr using reg_moments_mkt2, sort
drop _merge
sort mom_nbr
merge mom_nbr using reg_moments_mkt3, sort
drop _merge
sort mom_nbr
merge mom_nbr using reg_moments_mkt4, sort

drop mom_nbr _merge

outsheet using "reg_moments.csv", comma replace

*- Bootstrap and CIs
clear all
set more off, permanently

*---Mkt1

u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==1

cap program drop prog_regboot
program prog_regboot, eclass
tempname ests
	mat `ests'=J(1,21,.)

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==1 & couple==1, robust
mat `ests'[1,1]=_b[_cons]
mat `ests'[1,2]=_b[maleage]
mat `ests'[1,3]=_b[sqrage]
mat `ests'[1,4]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==2 & couple==1, robust
mat `ests'[1,5]=_b[_cons]
mat `ests'[1,6]=_b[maleage]
mat `ests'[1,7]=_b[sqrage]
mat `ests'[1,8]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==3 & couple==1, robust
mat `ests'[1,9]==_b[_cons]
mat `ests'[1,10]=_b[maleage]
mat `ests'[1,11]=_b[sqrage]
mat `ests'[1,12]=_b[K]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==1 & couple==1, robust
mat `ests'[1,13]=_b[_cons]
mat `ests'[1,14]=_b[X]
mat `ests'[1,15]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==2 & couple==1, robust
mat `ests'[1,16]=_b[_cons]
mat `ests'[1,17]=_b[X]
mat `ests'[1,18]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==3 & couple==1, robust
mat `ests'[1,19]=_b[_cons]
mat `ests'[1,20]=_b[X]
mat `ests'[1,21]=_b[sqr_X]

	mat list `ests'
	ereturn post `ests'
end

prog_regboot
bootstrap _b, reps(1000) seed(1218) saving(regboot_mat_1,replace): prog_regboot
clear
u regboot_mat_1, clear
forvalues j=1/21{
rename _b_c`j' regcoeff_`j'_1
}
gen rep=_n
sort rep

save reg_boot_1.dta, replace

*---Mkt2
u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==2

cap program drop prog_regboot
program prog_regboot, eclass

tempname ests
	mat `ests'=J(1,21,.) 

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==1 & couple==1, robust
mat `ests'[1,1]=_b[_cons]
mat `ests'[1,2]=_b[maleage]
mat `ests'[1,3]=_b[sqrage]
mat `ests'[1,4]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==2 & couple==1, robust
mat `ests'[1,5]=_b[_cons]
mat `ests'[1,6]=_b[maleage]
mat `ests'[1,7]=_b[sqrage]
mat `ests'[1,8]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==3 & couple==1, robust
mat `ests'[1,9]==_b[_cons]
mat `ests'[1,10]=_b[maleage]
mat `ests'[1,11]=_b[sqrage]
mat `ests'[1,12]=_b[K]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==1 & couple==1, robust
mat `ests'[1,13]=_b[_cons]
mat `ests'[1,14]=_b[X]
mat `ests'[1,15]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==2 & couple==1, robust
mat `ests'[1,16]=_b[_cons]
mat `ests'[1,17]=_b[X]
mat `ests'[1,18]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==3 & couple==1, robust
mat `ests'[1,19]=_b[_cons]
mat `ests'[1,20]=_b[X]
mat `ests'[1,21]=_b[sqr_X]

	mat list `ests'
	ereturn post `ests'
end

prog_regboot
bootstrap _b, reps(1000) seed(1218) saving(regboot_mat_2,replace): prog_regboot
clear
u regboot_mat_2, clear
forvalues j=1/21{
rename _b_c`j' regcoeff_`j'_2
}
gen rep=_n
sort rep

save reg_boot_2.dta, replace

*---Mkt3

u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==3

cap program drop prog_regboot
program prog_regboot, eclass

tempname ests
	mat `ests'=J(1,21,.) 

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==1 & couple==1, robust
mat `ests'[1,1]=_b[_cons]
mat `ests'[1,2]=_b[maleage]
mat `ests'[1,3]=_b[sqrage]
mat `ests'[1,4]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==2 & couple==1, robust
mat `ests'[1,5]=_b[_cons]
mat `ests'[1,6]=_b[maleage]
mat `ests'[1,7]=_b[sqrage]
mat `ests'[1,8]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==3 & couple==1, robust
mat `ests'[1,9]==_b[_cons]
mat `ests'[1,10]=_b[maleage]
mat `ests'[1,11]=_b[sqrage]
mat `ests'[1,12]=_b[K]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==1 & couple==1, robust
mat `ests'[1,13]=_b[_cons]
mat `ests'[1,14]=_b[X]
mat `ests'[1,15]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==2 & couple==1, robust
mat `ests'[1,16]=_b[_cons]
mat `ests'[1,17]=_b[X]
mat `ests'[1,18]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==3 & couple==1, robust
mat `ests'[1,19]=_b[_cons]
mat `ests'[1,20]=_b[X]
mat `ests'[1,21]=_b[sqr_X]

	mat list `ests'
	ereturn post `ests'
end

prog_regboot
bootstrap _b, reps(1000) seed(1218) saving(regboot_mat_3,replace): prog_regboot
clear
u regboot_mat_3, clear
forvalues j=1/21{
rename _b_c`j' regcoeff_`j'_3
}
gen rep=_n
sort rep

save reg_boot_3.dta, replace

*---Mkt4
u data_for_earnings-moments-allmkts.dta, clear
keep if USCBregion==4

cap program drop prog_regboot
program prog_regboot, eclass
tempname ests
	mat `ests'=J(1,21,.) 

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==1 & couple==1, robust
mat `ests'[1,1]=_b[_cons]
mat `ests'[1,2]=_b[maleage]
mat `ests'[1,3]=_b[sqrage]
mat `ests'[1,4]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==2 & couple==1, robust
mat `ests'[1,5]=_b[_cons]
mat `ests'[1,6]=_b[maleage]
mat `ests'[1,7]=_b[sqrage]
mat `ests'[1,8]=_b[K]

qui reg l_ly K maleage sqrage if man!=.  & maleduc_cat==3 & couple==1, robust
mat `ests'[1,9]==_b[_cons]
mat `ests'[1,10]=_b[maleage]
mat `ests'[1,11]=_b[sqrage]
mat `ests'[1,12]=_b[K]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==1 & couple==1, robust
mat `ests'[1,13]=_b[_cons]
mat `ests'[1,14]=_b[X]
mat `ests'[1,15]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==2 & couple==1, robust
mat `ests'[1,16]=_b[_cons]
mat `ests'[1,17]=_b[X]
mat `ests'[1,18]=_b[sqr_X]

reg lfem_ly X sqr_X if woman!=.  & femeduc_cat==3 & couple==1, robust
mat `ests'[1,19]=_b[_cons]
mat `ests'[1,20]=_b[X]
mat `ests'[1,21]=_b[sqr_X]

	mat list `ests'
	ereturn post `ests'
end

prog_regboot
bootstrap _b, reps(1000) seed(1218) saving(regboot_mat_4,replace): prog_regboot
clear
u regboot_mat_4, clear
forvalues j=1/21{
rename _b_c`j' regcoeff_`j'_4
}
gen rep=_n
sort rep
save reg_boot_4.dta, replace

*** Append all markets together

use reg_boot_1.dta, clear
merge rep using reg_boot_2
drop _merge
merge rep using reg_boot_3, sort
drop _merge
merge rep using reg_boot_4, sort
drop rep
drop _merge

outsheet using "reg_boot.csv", comma replace
}

*%--------------- Untargeted bootstrapped data moments and CIs ----------------%

clear all
set more off, permanently
u data_for_untargeted_moments, clear								  

qui{
*- Duration							  
bysort couple_type: sum duration if couple_marries==1 &ever_divorced==1
forvalue g=1(1)9{
gen duration_`g' =.
sum duration if couple_marries==1 &ever_divorced==1 & couple_type==`g'
replace duration_`g' = r(mean) if couple_marries==1 &ever_divorced==1 & couple_type==`g'
}

*- Divorce 
forvalues i=1(1)4{
gen div_hazard_I`i'=.
qui sum dh_I`i' if couple_type>=1 & couple_type<=3 
replace div_hazard_I`i'=r(mean) if couple_type>=1 & couple_type<=3
qui sum dh_I`i' if couple_type>=4 & couple_type<=6 
replace div_hazard_I`i'=r(mean) if couple_type>=4 & couple_type<=6
qui sum dh_I`i' if couple_type>=7 & couple_type<=9 
replace div_hazard_I`i'=r(mean) if couple_type>=7 & couple_type<=9
}

*- Stay-at-home wife 
forvalues i=1(1)4{
gen hw_I`i'=.
qui sum sahw_I`i' if couple_type>=1 & couple_type<=3 
replace hw_I`i'=r(mean) if couple_type>=1 & couple_type<=3
qui sum sahw_I`i' if couple_type>=4 & couple_type<=6
replace hw_I`i'=r(mean) if couple_type>=4 & couple_type<=6
qui sum sahw_I`i' if couple_type>=7 & couple_type<=9 
replace hw_I`i'=r(mean) if couple_type>=7 & couple_type<=9
}
forvalues i=1(1)4{
forvalues g=1(1)9{
gen hw_I`i'_`g'=.
qui sum sahw_I`i' if couple_type==`g' 
replace hw_I`i'_`g'=r(mean) if couple_type==`g' 
}
}

forvalues i=1(1)10{
gen hw_t`i'=.
qui sum sahw_t`i' if couple_type>=1 & couple_type<=3 
replace hw_t`i'=r(mean) if couple_type>=1 & couple_type<=3
qui sum sahw_t`i' if couple_type>=4 & couple_type<=6 
replace hw_t`i'=r(mean) if couple_type>=4 & couple_type<=6
qui sum sahw_t`i' if couple_type>=7 & couple_type<=9 
replace hw_t`i'=r(mean) if couple_type>=7 & couple_type<=9
} 
								  
* Outsheet data
keep couple_type duration_1-duration_9 hw_I1_1-hw_I4_9
drop if couple_type>10
egen duration = rowtotal(duration_1-duration_9)
drop duration_1-duration_9
drop if duration==0
egen sahw_1= rowtotal(hw_I1_1-hw_I1_9)
egen sahw_2= rowtotal(hw_I2_1-hw_I2_9)
egen sahw_3= rowtotal(hw_I3_1-hw_I3_9)
egen sahw_4= rowtotal(hw_I4_1-hw_I4_9)
keep couple_type duration - sahw_4
duplicates  drop
sort couple_type
outsheet  using "data_untargeted.csv", comma replace
								  
* Bootstrap
clear all
set more off, permanently
u data_for_untargeted_moments								  

cap program drop prog_untg_momentboot
program prog_untg_momentboot, eclass
tempname ests
	mat `ests'=J(1,9+9+9+9+9,.) 
*- Duration
forvalues g=1(1)9{
qui sum duration if couple_marries==1 &ever_divorced==1 & couple_type==`g'
mat `ests'[1,`g']=r(mean)
}
*- Stay at home period 1
forvalues g=1(1)9{
qui sum sahw_I1 if couple_type==`g'
mat `ests'[1,`g'+9]=r(mean)
}
*- Stay at home period 2
forvalues g=1(1)9{
qui sum sahw_I2 if couple_type==`g'
mat `ests'[1,`g'+9 +9]=r(mean)
}
*- Stay at home period 3
forvalues g=1(1)9{
qui sum sahw_I3 if couple_type==`g'
mat `ests'[1,`g'+9 +9 +9]=r(mean)
}

*- Stay at home period 4
forvalues g=1(1)9{
qui sum sahw_I4 if couple_type==`g'
mat `ests'[1,`g'+9 +9 +9 +9]=r(mean)
}
	ereturn post `ests'
end

prog_untg_momentboot
bootstrap _b, reps(1000) seed(1218) saving(untg_momentboot_mat,replace): prog_untg_momentboot

*- Create 95 CI bounds and export
clear
use untg_momentboot_mat
collapse (p5) _b_c1- _b_c45

save untg_moments_boot_lb.dta, replace

clear
use untg_momentboot_mat
collapse (p95) _b_c1- _b_c45
append using untg_moments_boot_lb.dta

save untg_moments_boot_CI.dta, replace

outsheet using "untg_moments_boot_CI.csv", comma replace

}



























